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Definite integrals of analytic functions can often be evaluated efficiently 
by the trapezoidal rule after a suitable transformation. Here the work of 
Moran 1 and Schwartz 2 along this line is extended. First the dependence of 
the error on the spacing is discussed, and then several types of transforma- 
tions are described and applied to integrals of technical interest. 

I. INTRODUCTION 

Quite often the problem of determining the value of a definite in- 
tegral arises. When the integral cannot be readily evaluated by analysis, 
we must resort to numerical methods. Here we discuss a method of 
numerical quadrature which gives promise of being useful in evaluating 
some types of integrals that are difficult to handle by conventional 
numerical methods. 

In particular, we consider the problem (Moran 1 and Schwartz 2 ) of 
transforming a given integral of an analytic function f{x) into a rapidly 
converging one (with limits ± °c ) which can be efficiently evaluated by 
the trapezoidal rule, 

[ m f(x)dx = h t f(nh) - E. (1) 

J— oo n=— oo 

The integral and series are assumed to converge. In addition to the 
trapezoidal error E, a second error is introduced when the series is 
truncated in the process of computation. It is supposed that both 
errors are made negligible, E by taking h small, and the truncation 
error by taking enough terms in the series. The feature which makes the 
use of (1) attractive is that E often decreases in proportion to 
exp ( — C/h) as h decreases, C being a constant. Thus if h gives three- 
figure accuracy, h/2 will give six-figure accuracy in many cases. 
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The transformations, i.e., the changes of variable of integration used 
to carry the limits of the given integral into ± w , are usually con- 
structed by combining functions which are readily computed, such as 
powers and exponential functions. 

Schwartz 2 recommends the following procedure for evaluating the 
integral of an analytic function f(u) : (i) change the variable from u to 
v so as to make the integration with respect to v extend from — °° to 
+ oo t (u) make the further transformation 

v = e 1 — e~ x 

to increase the rate of convergence, and then (Hi) evaluate the trape- 
zoidal sum, truncating when contributions fall below the desired ac- 
curacy, and reduce the spacing h until the answer has the desired 
accuracy. 

Here we present a summary of several variations of Schwartz's pro- 
cedure. Details are given in a report by the author. 3 The dependence of 
the error E on the spacing h is first reviewed, and then examples are 
used to illustrate the evaluation of various types of integrals. 

II. THE DEPENDENCE OF E ON fl 

The trapezoidal error E can be expressed in several ways. For ex- 
ample, it is the remainder in the Euler-Maclaurin sum formula. Again, 
it can be written as the sum of contour integrals with integrands 
/(z)/[exp {±i2tz/K) - 1] (Ref. 4, p. 145, Problem 7, and Refs. 5, 6, 
7, 8, and 9). Here we follow Fettis 5 and use Poisson's summation for- 
mula which, when applied to (1), gives 

E = ( £ +f)f /(z) exp (i2irzk/h)dz. (2) 

Let f(z) be analytic throughout a strip in the z-plane containing the 
real z-axis and assume that suitable convergence conditions are satis- 
fied. Then the paths of integration in the terms for k > in (2) can be 
displaced upwards to make Im (z) > 0. It follows that | exp (i2irkz/h) \ 
= exp ( — 2-rrk Im (z)/h) becomes small when h becomes small and z is 
on the path. Furthermore, as h — » the terms for k > 1 become 
negligible in comparison with the term for k = 1. A similar argument 
holds for the k < terms, and as h — > we have the asymptotic result 

E ~ R+ + R- 

where R+ and R- are the k = 1 and k = - 1 terms, respectively, in 
(2). For the important case in which /(z) is real on the real z-axis, R- 
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is equal to the conjugate complex R' + of R + and E is given asymptoti- 
cally by 

E - R + + RT+, 



R 



;+ = \™ f(z) exp {i2irz/h)dz. (3) 

When h is small, the trapezoidal error E given by (3) can be viewed 
as the sum of contributions from singularities of f(z) and saddle points 
of f(z) exp (i2wz/h). At the saddle points the derivative d<p(z)/dz is 
zero, <p(z) being defined to within a multiple of 2rt by exp [*>(«)] 
= /(z) exp (i2irz/h). 

This picture of £ is suggested by the following remarks. The path of 
integration in (3) can be deformed upwards towards z = too in the 
complex z-plane until it becomes an optimal path comprised of paths 
of steepest descent and ascent passing through one or more saddle 
points. The path runs from — <*; to + w and may have detours running 
out to, and returning from, infinity. It may also have loops around 
some of the singularities of /(z). When h is small, the factor exp (i2*z/h) 
decreases rapidly as Im (z) increases, and the only significant contribu- 
tions to R+ come from the portions of the path near the singularities 
and near the saddle points not associated with singularities. 

An approximate expression for a typical contribution can be obtained 
by expanding the contribution about the corresponding singularity or 
saddle point and taking the leading term. Such expansions are usually 
asymptotic in h. For estimating orders of magnitude we can use the 
dominant factors in the leading terms: 

(Contribution to R+ of a saddle point at z ) « exp C<^(z )], 
(Contribution to R+ of a singularity at Z\) ~ exp (i2wZi/h). 

As h decreases, E may either decrease steadily or may oscillate with 
decreasing amplitude depending upon how the dominant contributions 
combine. 

If there are no singularities or saddle points in the finite part of the 
z-plane, the trapezoidal rule may give the exact value of the integral 
when h is less than some fixed value. This is associated with the 
sampling theorem for band-limited functions. For example, if m and n 
are positive integers such that m — n = or 2, 4, • • • , the integral 

/ = / sin m xdx/x n (4) 

is exactly equal to the trapezoidal sum when h < 2ir/m (h can equal 
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2tt/w if n 2i 2). The proof follows from the fact that all of the terms 
in the series (2) for E vanish when h < 2tr/m, as can be shown by de- 
forming the paths of integration into infinite semicircles and using 
Jordan's lemma (Ref. 4, p. 115). As a check, note that for m = n = 1 
or 2 we can take h = v. Then I = ir because there is only one nonzero 
term in the sum. 

The foregoing discussion shows that the structure of E can be de- 
termined by computing the saddle points and associated paths of 
steepest descent for /(z) exp (i2irz/h). This is done in the report 3 for 
the examples (G) and (25) given below. The path for example (6) is 
shown in Fig. 1 and discussed in the appendix. However, computations 
of this sort are laborious. In practice it appears that the dependence of 
E on h is most easily determined by computing the trapezoidal sum for 
a sequence of decreasing values of h, bearing in mind the possibility 
that E may go through zero for some values of h. 

Incidentally, arguments similar to those given in this section show 
that the trapezoidal rule also works well when it is used to evaluate in- 
tegrals of periodic analytic functions in which the integration extends 
over a period. 

III. CONTRIBUTIONS TO E — EXAMPLE 

Goodwin 6 has pointed out that the trapezoidal rule usually performs 
well for integrals of the type 



'-/-" 



g(x) exp ( — bx 2 )dx. 



(5) 




t5 



^5 



3 



--3 



z- PLANE 




Fig. 1 — Steepest descent paths for exp O(z)] when a = 2.4 and h = 0.8 in both 
(30) and example (6). The points z , z is , zi> are saddle points; z,, z 2 are branch points. 
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Computations 3 made with 6=1, 

g(x) = (x* + a 4 )-*, (6) 

and a = 2.4 show that, as h decreases, E first decreases steadily and 
then near h = 1 (where E/I « 10 -4 ) E starts to oscillate with decreas- 
ing amplitude. This behavior can be explained in terms of a contribu- 
tion to E of approximately 

2tt4 Re lgiiir/h) exp (-**/*»)] (7) 

from a saddle point near z = iir/h (see Fig. 1, eq. (32), and Goodwin 8 ) 
and contributions from the branch points at Z\ = a exp (zV/4) and 
zi = a exp (2'37r/4). The expressions for the branch point contributions 
are somewhat more complicated than (7), but all that need be noted 
here is that they contain exp {i2irz k /h), k = 1, 2, as a dominant factor. 
Adding the three contributions and neglecting multipliers such as 
g(iw/h) in (7) shows that E is roughly 

exp (-ir 2 /h 2 ) + (cos a) exp (-0) (8) 

where /3 = 2-'ira/h and cos a oscillates with increasing frequency as h 
decreases. The steady decrease of E with h, dominated by exp ( — Tr 2 /h 2 ), 
changes to an oscillating decrease when /3 = w 2 /h 2 . Solving for h with 
a = 2.4 gives h = 0.93, which agrees with the observed h « 1.0. 

When g(x) in (5) is algebraic and g(z) has no singularities inside the 
rectangle with corners at ±z , ±|2o|, where z = ir/(bh) and Re (6) 
> 0, the error E tends to be dominated by the saddle point contribu- 
tion. This contribution is approximately | exp (bzl) | in the sense that 
exp ( — ir 2 /h 2 ) approximates the approximation (7). When the rectangle 
contains singularities, their contributions dominate. For the example 
(6), zo = iv/h and the rectangle becomes a square which expands as h 
decreases. The behavior of E changes when the sides of the square 
sweep across the branch points. 

IV. INTEGRALS WITH BOTH LIMITS FINITE 

Consider the integral 

I = f (u - a) a -'(6 - uy- l f(u)du, a, 13 > 0, (9) 

where f(u) is analytic and 0(1), f(a) and f(b) =^ 0, and a and b are 
finite. The transformation 

u = (be'' + ae">)/(e v + e-), . 

du/dv = 2(6 - a)/(e v + e-') 2 , ( ; 
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carries the limits into v = =fc «> . The associated equations 

u- a = e v (b - a)/(e v + e~ v ), 
b - u = e~ v (b - a)/(e" + e~ v ), 

and du/dv show that the dominant factors in the integrand when 
v — > -f oo and v — * - °° are exp ( — 2@v) and exp (2av), respectively. 
One might expect, and it is confirmed by computation, that in a good 
transformation the final integral should converge at — °° at nearly the 
same rate as it does at + <» . Therefore the further change of variable 

v = c(P~ l e x — a -1 e _x ) 
dv/dx = c(/3-'e x + or l <T*) 

is made to equalize the rates of convergence at x = ± °o . The transfor- 
mation (11) makes the integrand behave roughly as exp [-2c exp |x|] 
as x — > ± °° when the effect of f(u) is ignored. 

The constant c in (11) can be chosen somewhat arbitrarily. It is 
helpful, but not necessary, to have 

c ^ 7r(/3a)V4. (12) 

The inequality (12) for c guarantees that the singularities in the com- 
plex x-plane due to the vanishing of e" + e~" are at least tt/2 distant 
from the real z-axis. It says nothing about the singularities and saddle 
points introduced by f(u). 

Thus the integral to be evaluated by the trapezoidal rule is 

1 - L {u " a)a ~ Kb " uy ~ lJ{u) Tv i dx (13) 

which can also be written as 

/oo p(a—0)v /7ji 

Here x is the variable of integration and, in writing the program, u, v, 
du/dv, and dv/dx are given by (10) and (11). 

As an example of (9) and (13), consider the beta function 

/ = J [sin w]"- 1 sin ( ^ - uj\ du 
= £ [sin I.]-' [sin ( | - i,)]*" 1 J | d* 
= «/_" [sin y]-« [sin (| - «)]'" | L*+JL^ ] dx 
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where the third line is to be evaluated by the trapezoidal rule. For 
a = 0.95 and /3 = 0.05, the value of / is known to be 20.748 732 • • • 
and the inequality (12) for the multiplier c gives c ^ 0.171. 

Computations show that the values 0.171, 0.1, and 0.05 for c all 
give six-figure accuracy with h = 0.5 and about 20 terms in the trape- 
zoidal sum (e.g., for c = 0. 1, h = 0.5, and 21 terms, the computed 
value is 20.748 729). When c = 1, about 70 terms (with h = 0.075) are 
required to achieve the same accuracy. 

Instead of computing the a — 1 power of simply [sin u~\, it was 
found better to compute the a — 1 power of [ie v + e-') 2 sin u~\, and 

similarly for sin ( ~ — u J • In general, it is usually helpful to combine 

as much as possible of the l/(e" + e -1 ') 2 contained in du/dv with other 
factors in the integrand in order to avoid underflow and overflow. 

When a = /3 = 1, the transformations (10) and (11) reduce to the 
ones used by Schwartz 2 except for the coefficient c ^ x/4 = 0.785. To 
illustrate this case take (Kajfez 1 ") 

/ = (-tt/40) f ° exp (u/4) sin [0.4tt exp (u/4)ldu 

JlQ 

in which the integrand oscillates through about 6 cycles. Using (10) 
and (11) with a = 10, 6 = 15, c = 0.785 shows that the trapezoidal 
rule with h = 0.09 and 60 terms gives / = —0.0195495 compared with 
the true value —0.0195488 • • • . For relative errors less than about 0.01, 
the trapezoidal rule requires less terms than the spline quadrature 
methods considered in Ref. 10, but this is offset somewhat by the more 
complicated terms introduced by (10) and (11). 

V. INTEGRALS WITH LIMITS 0, °o CONTAINING W a_1 (l + u)~ a ~ P 

The integral 

I = T M o-i(i + u)-°-ef(u)du, a, > 0, (14) 

Jo 

where /(0) 9^ and f(u) is analytic and 0(1), can be handled by the 
transformations 

u = e v , v = c(@- l e x — a-^-') (15) 

where c ^ 7r(a/3)V2. For a = 3, /3 = 2 the inequality for c becomes 
c ^ 3.85, and for a = 0.2, = 0.1 it becomes c ^ 0.222. Computations 
were made for these values of a and (3 with f(u) = 1. Values of h and 
the number of terms N in the trapezoidal sum required for seven- 
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figure accuracy were 


found to b 


e as foil 


3\vs: 




a 





c 


h 


N 


3.0 


2.0 


3.85 


0.25 


15 






2.00 


0.35 


15 






5.00 


0.10 


40 


0.2 


0.1 


0.22 


0.45 


25 






0.08 


0.45 


25 






0.45 


0.25 


35 



VI. INTEGRALS WITH LIMITS 0, °° CONTAINING U a l exp (~u) 

For the integrals 

I = f u°- 1 e-' t f(u)du, a > 0, (16) 

Jo 

where /(0) ^ and f(u) is analytic and 0(1), we can use 

u = e v , v = x — a~ l e~ x . (17) 

Computations for the case f(u) = 1 and a = 1 gave the following 

results: 

h N = No. of terms Trap, values of / 

0.4 15 0.9999 9997 

0.0 10 0.9999 8711 

0.8 7 0.9998 2442 

Repeating the computations with u = c exp [x — c exp (—x)] and 
c = 0.5, 2.0, and 4.0 showed that the magnitude of the error depends 
only slightly on c. 

VII. INTEGRANDS WHICH CHANGE RAPIDLY NEAR A POINT 

When the integrand contains a factor, say F(t) where t is the variable 
of integration, which changes rapidly near a point it is sometimes help- 
ful to change to a new variable of integration u where du/dt = F(t) 
and the constant of integration is chosen at our convenience. The suc- 
cess of the transformation depends upon the ease of inverting to get 
t as an easily computed function of u. 

As an example, consider 

I = J' e'(t 2 + a*)-Wt (18) 

where a is small (Smith and Lyness 11 )- Taking du/dt = F(t) 
= (t 2 + a 2 )~ { , w = arcsinh (t/a), and t = a sinh u carries (18) into an 
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integral of the form (9) with b = —a = A and a = = 1 : 

du dv 
■ , dv dx 



I = f e'du = I e l — -' — dx 



= AAc ( X e'(e' + e~*)dx/(e v + e~") 2 . 

J — 00 

Here A = arcsinh (1/a), t = a sinh (A tanh v), v = 2c sinh x, and 
c ^ 7r/4 = 0.785. Computations with c = 0.3, A = 0.2, and 40 terms 
in the trapezoidal sum show that / = 29.538 618 ± 10 -6 when 
a = 10- 6 . 

If the rapidly changing factor has the more general form 

F(t) = (t> + a 2 )- 9 , 

the transformation t = a sinh u can still be used to carry the integral 
into the form (9), but computation shows that the trapezoidal rule re- 
quires more terms as d moves away from \. For example, when the 
exponent — \ in (18) is replaced by — f, i.e., 6 = f, computations with 
a = 10~ 6 , c = 0.785, and h = 0.03 show that 100 terms give the value 
5240.808 for / whereas the true value is 5240.806 

VIII. THE FERMI-DIHAC INTEGRAL 

The Fermi-Dirac integral (tabulated by Blakemore 12 ) 

1 = rh)J* ta ~ ldt/(l + e '~ a) > a>0 > (19) 

has an integrand which changes rapidly near t = a when a is large. 
Section VII suggests taking du/dt = F(t) = 1/(1 + e' - °). Choosing 
the constant of integration to make u = at t = gives 

u = tn(\ + e~ a ) - fn(e~ l + e~ a ), 

t = -fn(e~" + e-"~ a - e~ a ), (20) 



where b = (n(l + e a ). When u tends to 0, t — > (I -J- e~ a )u, and hence 
the integi-al (20) is of the form (9) with /3 = 1. 
The transformations (10) and (11) carry / into 

, 1 f* . .du dv , 

I (a) y-oo dvdx 

= W^r ta ~ l ( e ' + a-'e-^dx/ie" + <r") 2 (21) 

I (a) J- K 
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where c ^ 0.7850*, t is defined in terms of u by (20), and 

u = be v /(e" + e~ v ), v = c(e r — oT l «T x ). 

The integral (21), with c = 0.5, was used to compute / for a = \ 
and f with a between and 20. Difficulty in computing t for small 
values of u was avoided by using three terms in the expansion of 
ln{\ - z) when z = exp (u — a) — exp (-a) was less than 0.001. The 
following tabulation shows the results for the typical values a = \ and 
a = 10. 

h N = No. of terms Trap, values of I 

0.2 34 3.5527 792 

0.3 22 792 

0.4 17 795 

0.5 14 742 

The above transformation has been used by W. K. Kent in a study 
of the charge distribution in a charge coupled device and I am indebted 
to him for helpful discussions regarding his experience with (21). 

IX. INTEGRALS WITH LIMITS 0, <*> AND OSCILLATING INTEGRANDS 

Kluyver's Bessel function (random walk) integral for the probability 
P that the resultant of the sum of m randomly phased unit vectors in 
a plane be less than r in length is (Bennett 13 , Greenwood and Durand 14 ) 

P = f a r Ji (rw)[Jo (u)2 m du. (22) 

Jo 

This integral is typical of a class of integrals that are difficult to evalu- 
ate by any means. They are characterized by rather slow convergence 
and an integrand which tends to oscillate at a regular rate as u -* « . 
In this section we consider the evaluation of such integrals by the 
trapezoidal rule when the rate of convergence is not too slow. 

Some general remarks can be made concerning integrals that behave 
like (22). In order that 

I = T f(u)du 
Jo 

may represent the typical integral of this section, f(u) must tend, as 
u — >«>, to a form that can be written as a steadily decreasing factor 
times the sum of a finite number of sinusoidal terms whose periods tend 
to constant values (as u— >•»). Define h to be the shortest con- 
stant period (so that the most rapidly oscillating term varies as 
cos [(2ttm/A ) + 0]). The interval h is related to the sampling theorem 



INTEGRALS OF ANALYTIC FUNCTIONS 717 

for band-limited functions in that l/h plays the role of the "band- 
width" of f(u) at u = oo . Quite often the required accuracy in / can 
be obtained by using values of h which are close to h . 

If the integrand is an even function of u, the trapezoidal rule can be 
applied directly. If the integrand is not even, the integral can be evalu- 
ated by setting 

u = a(n{\ + e' la ), du/dx = e x ' a /{\ + e zla ) (23) 

and then using the trapezoidal rule. Computations described below 
show that the choice a = 1 works well for (22). More will be said later 
about the choice of a. The transformation (23) takes advantage of the 
fact that the trapezoidal error E is often small, or zero, for regularly 
oscillating integrands. Although many terms may be needed in some 
cases, the present method compares favorably with competing ones. 

The choice of a in (23) depends upon the behavior of f(u) near u = 0. 
Suppose that f(n) tends to Cu" as u — > 0. Here C is a constant and 
v > —\. After the change of variable (23), the integrand is a function 
of x which approaches as exp [_{v -+- l)x/a] when x — > — °o . When 
(v -f- \)h/a is too small, successive terms in the negative x portion of 
the trapezoidal sum are nearly equal and an unduly large number of 
terms must be taken to achieve a small truncation error on the left. 
When (v -f- l)h/a is too large, the successive terms differ by a large 
amount and the trapezoidal error E tends to be large. The problem is 
to choose a value of a which balances these two effects and at the same 
time allows h to be large. The choice a = (v + l)/i works well for all of 
the cases that have been tried. 

Some insight regarding a good choice of h for (22) can be obtained as 
follows. Consider the integral, say K, obtained by replacing Ji(ru) by 
Jo(ru) in (22) and taking the limits of integration to be u = ± °o . With 
the help of the asymptotic expression for Jo(z) and the procedure used 
to deal with (4), it can be shown that when K is evaluated by the trape- 
zoidal rule the error is zero if h < h where h = 2w/(r + wi). Further- 
more, when P is transformed by (23) and then evaluated by the trape- 
zoidal rule the error E is relatively small when h is only slightly less 
than h — as might be hoped from the similar behavior of the integrands 
in P and K as u — ><» . A saddle point analysis of R+ in (3) leads to the 
rough estimate 

\E\ « 2r*(2ir 2 )~ (w+1)/2 exp [>(r + m - 2irhr l )2 (24) 

for the error in the trapezoidal sum for P when a = 1 . 

In the example (22) the asymptotic expression for the integrand con- 
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tains the product cos [ru — (3tt/4)] cos ra [u - (tt/4)] which can be 
written as the sum of m + 1 sinusoidal terms. The most rapidly oscil- 
lating term is proportional to cos \_(r + m)u — (m + 3) (tt/4)] and 
the quantity h Q = 2w/(r + m) now appears as the shortest period. 

Now we turn to the details of the evaluation of the integral (22) for 
P. Let r and m have the representative values r = 4 and m = 6. Then 
h = 2w/(r + w) = 0.628. Since J (z) -» 1 and Ji{z) -> 2/2 as z -> 0, 
the exponent c is 1. The suggested value a = (v + l)/io gives a = 2/i 
= 1.26, and since the choice of a is not critical, we take a = 1. Putting 
a = 1 in (23), substituting in (22), and using the trapezoidal rule gives 
the values of P shown in column 3 : 



h 


No. of terms 


P 


Error: Col. 3 


| E | from (24) 


0.475 


286 


0.9375 5485 






0.500 


272 


5475 


-10. X 10~ 8 


3.5 X 10~ 8 


0.525 


259 


5437 


- 4.8 X 10~ 7 


2.8 X 10- 7 


0.550 


247 


5354 


- 1.3 X 10-" 


1.3 X 10" 8 


0.575 


236 


5791 


+ 3.1 X 10- 6 


6.9 X 10- 6 


0.600 


226 


0.9375 9798 


+ 4.3 X 10"« 


2.4 X 10- 5 


0.625 


217 


0.9376 9974 


+ 1.4 X 10" 4 


0.94 X 10- 4 



The fourth column gives the error estimated from column 3. The 
fifth column shows the approximation (24) for \E\. The trapezoidal 
sum was truncated at x = 124. Beyond 124 the absolute value of the 
integrand remains less than 2 X 10 -8 and its amplitude decreases as 
r -7/2 ]sj ote that the error starts to be appreciable as h approaches the 
critical value h = 2ir/(r + m) = 0.628. 

X. THE INTEGRAL OF U k exp ( — M 2 — aW 1 ) FROM TO « 

The integral 

h(a) = ( K u*exp (-w 2 - au- l )du, Re (a) ^ 0, (25) 

Jo 

is of interest in some physical problems (I wish to thank J. N. Lyness 
for calling my attention to this example). First let A; be and a be posi- 
tive real. We seek a change of variable from u to x, with new limits 
x = ± oo, which will make u 2 tend to exp (2x) as x — »<» and a/'u tend 
to exp ( — 2x) as x — * — °o . This leads to 

u = ae*/(a + e~ x ). (26) 
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Substituting (26) in (25), taking the special case k = 0, a = 1, and 
applying the trapezoidal rule gives the values of 7 (1) shown in column 
3: 



h 


No. of terms 


UD 


Obs. error 


\E\ from (27) 


0.1 


38 


0.1500 4597 


2 X 10~ 8 




0.2 


19 


0.1500 4597 


2 X 10~ 8 


5 X 10~ 9 


0.3 


12 


0.1500 4835 


2.40 X 10~ 6 


7 X 10-« 


0.4 


10 


0.1501 2711 


S.12 X 10- 6 


20 X 10- 5 



The last column lists a rough approximation obtained by saddle point 
analysis : 

w-»[-i! + (t)T ™ 

The observed error shown in column 4 is the trapezoidal sum (column 
3) minus the value 0.1500 4595 of / a (l) computed from 

J ^-ia/^^ i)r ( il 4 ± - 1 )-* 

f (~ a )" r ( k-n+l \ f, t (-l)»+*o 2 "+* 
~ £„ 2(n!) V 2 ) + i,2" 



^2(2n + A;)!r(n4-i) 



+ £, n!(an + * + l)l ( j 

where c < min (0, A- + 1), \f/(x) = (dr(.r)/rf.r)/r(.T), and the series 
holds for k = — 1, 0, 1, 2, 3, ■ • ■ except that the first sum is omitted 
when k = — 1 . 

When a is complex, say a = pe' a where |a| ^ ir/2, we tilt the path 
of integration in (25) by setting u = ve' e where | d\ ^ ir/4 and |<x — 6\ 
< 7r/2. If we choose 6 to be a/3, the new integral 

h( P e ia ) = exp [i(* + l)«/3] f° v k exp[-(i> 2 + p»- I )c i2a/3 ]dw (28) 

Jo 

can be evaluated by using the substitution (26) with u and a replaced 
by v and p, and then applying the trapezoidal rule. For the physically 
important case of k = 3 and imaginary a [Is(ip) and Ik{a), k = 1, 2, 3, 
are tabulated in NBS Handbook, 15 Section 27.5]], the error can be kept 
below 1 X 10 -6 for a = t'0.001 by using h = 0.04 and 95 terms in the 
trapezoidal sum. As a increases to z'10.0, h can increase to 0.08 and the 
required number of terms decrease to 40. 
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APPENDIX 

Examples of Paths of Steepest Descent for R+ 

For the example (6), the integral (3) for R+ becomes 

R + = r exp Mz)ldz, (29) 

J — 00 

v (z) = _ 2 2 _j_ Qrk-b - \tn(z A + a 4 ). (30) 

The saddle point equation d<p(z)/dz = <p'(z) = is of the 5th degree 
in z. Solving by Newton's rule or otherwise gives 5 saddle points, 3 of 
which are in the half-plane Im (z) > 0. They are shown as small circles 
in Fig. 1 for the case a = 2.4 and h = 0.8. One is on the imaginary 
z-axis at z - «4.14, and the other two (z u , zi.) are at ±1.648 + tl.629 
near the branch points (z h z 2 ) at o(=bl + *)/2* = 1.697(±1 + i). 

The path of steepest descent through the saddle point z (a path on 
which Im [>(z) - v>(z„)] = 0) was computed by: (i) evaluating the 
phase of the coefficient [-27r/^"(z )] i in the saddle point contribution 

/exp[^(z)]dz ~ [-27r/^"(zo)?exp [^(z«)] (31) 

to determine the direction of the path through z n ; (ii) selecting two 
starting points (one for each branch of the path) on opposite sides of, 
but close to, z ; and (in) applying 

z ( = z,_i + d e - lf d t = - | <p'(z ( ) | &/<p'(zi) 

to compute the path step by step, A being the step length. The other 
paths of steepest descent shown in Fig. 1 were computed in the same 
way. Paths of steepest descent through a saddle point may run down 
into a "lower" saddle point, i.e., one having a more negative Re <p(z), 
or may end at a point (possibly z = °o ) where Re <p(z) = — <* . 

Figure 1 shows that the path of integration for R+ can be deformed 
in a natural way into three portions: the loop around z 2 (which encloses 
the branch cut from z 2 ), the path from - oo + ivhr 1 through z to 
4- oo + i-irkr 1 , and lastly the loop around Zi. The path directions are 
denoted by arrows. 

As h increases (from 0.8), z in Fig. 1 moves downward towards the 
origin. Eventually h reaches a critical value at which the path of 
steepest descent from z runs directly into z u and z 2 „. For still larger 
values of h, the loops around Zi and z 2 lie above the path through z , 
and the deformed path of integration for R + consists only of the path 
running from - °° + zVA -1 to + « 4- iwhr 1 through z . The sudden 
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change in the path as h passes through the critical value is related to 
Stokes phenomena in the theory of asymptotic expansions. 

The approximation (7) for the contribution of z to E ~ 2 Re (R+) 
can be obtained by approximating the saddle point equation <p'{z) = 
by <Pa(z) = where 

<Pa (z) = -z 2 + i2irh~ l z (32) 

is the most important part of the expression (30) for <p(z) near z = z . 
The equation tp' A (z) = gives the approximation z A = iir/h for z . 
Substituting <p' A \z A ) = — 2 and <p(z A ) from (30) in place of ¥>"(z ) and 
<p(zo), respectively, in the saddle point contribution (31) to R + leads 
to (7). 

Figure 2 shows the paths of steepest descent and the path of integra- 
tion used to estimate E for the integral (25), 7 (1), when h = 0.2. For 
A: = and a = 1, the v(z) in the integral (29) for R + is 

+ 2) 



<p{z) = -u 1 - u~ l -V i2irzh- 1 + In 



\ e 2 '(e- + 2) ] 
L («' + I) 2 J 



(33) 



where u is given by (2G) with z in place of x. The deformed path of in- 
tegration for R + shown in Fig. 2 contains the arbitrary bridging seg- 
ment AB and shows that the saddle points z x and z a are the main con- 
tributors to R+. 




Fig. 2 — Steepest descent paths for exp [^(z)] when h = 0.2 in (33). The arrows 
mark a deformed path of integration for R + corresponding to 7<,(1) in (25). The points 
Z\, Zt, • • • , z g are saddle points. 
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Finally, if one wishes to integrate along a path of steepest descent, 
a combination of the above path computation and the trapezoidal 
rule with spacing A suggests itself. However, greater accuracy can be 
achieved by either (i) first computing the entire path, approximating 
it by several straight-line segments (sometimes one will do), and using 
Romberg integration on each segment, or (n) taking the step size A 
relatively large in the path computation, and then using Romberg 
integration over each linear segment of length A (by dividing d ( into 
lengths of, say, A/8). 
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